

# Here I redefine the getbetas function to get 
getbetas=function(object,row,column)
{
  if (any(class(object)=="eiMD")){
    Beta=object$draws$Beta
    if (is.null(Beta))
      stop("blablabla")
  }
  else if (any(class(object)=="eiMD.beta")){
    Beta=object
  }
  else stop("blablablabla")
  if (is.mcmc(Beta)){
    tnames=strsplit(colnames(Beta),"beta.")
    get2=function(x) x[2]
    idx=strsplit(sapply(tnames,get2),".",fixed=TRUE)
    idx=as.list(as.data.frame(matrix(unlist(idx),byrow=TRUE,
                                     nrow=length(idx),ncol=length(idx[[1]]))))
    idx=lapply(idx,as.character)
    idx=lapply(idx,unique)
    idx[[4]]=1:length(Beta[,1])
    names(idx)=c("rows","columns","precints","sims")
    Betas=array(t(object$draws$Beta),dim=sapply(idx,length),dimnames=idx)
  }
  else {
    idx=dimnames(Beta)
  }
  if (missing(row) | missing(column)) {
    stop("blablabla")
  }
  frow=idx$rows[grep(row,idx$rows)]
  fcolumn=idx$columns[grep(column,idx$columns)]
  if (length(frow)<1)
    stop("blablabla")
  if (length(frow)>2)
    stop("blablabla")
  if (length(fcolumn)<1)
    stop("blablabla")
  if (length(fcolumn)>2)
    stop("blablabla")
  usebetas=t(Betas[frow,fcolumn, ,])
  quant.low=function(x) {
    quantile(x,c(0.5-0.95/2))
  }   
  quant.high=function(x) {
    quantile(x,c(0.5+0.95/2))
  }       
  est=apply(usebetas,2,mean)
  low=apply(usebetas,2,quant.low)
  high=apply(usebetas,2,quant.high)
  output=data.frame(data$mesa,low,est,high)
  return(output)
}
